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ABSTRACT 

We show that current clustering observations of quasars and luminous AGN can be explained by a merger 
model augmented by feedback from outflows. Using numerical simulations large enough to study clustering 
out to 25 comoving Mpc, we calculate correlation functions, biases, and correlation lengths as a function 
of AGN redshift and optical and X-ray luminosity. At optical wavelengths, our results match a wide range 
of current observations and generate predictions for future data sets. We reproduce the weak luminosity de- 
pendence of clustering over the currently well-measured range, and predict a much stronger dependence at 
higher luminosities. The increase in the amplitude of binary quasar clustering observed in the Sloan Digital 
Sky Survey (SDSS) is also reproduced and is predicted to occur at higher redshift, an effect that is due to the 
one halo term in the correlation function. On the other hand, our results do not match the rapid evolution of the 
correlation length observed in the SDSS at z ~ 3, a discrepancy that is at least partially due to differences in the 
scales probed by our simulation versus this survey. In fact, we show that changing the distances sampled from 
our simulations can produce changes as large as 40% in the fitted correlation lengths. Finally, in the X-ray, our 
simulations produce correlation lengths similar to that observed in the Chandra Deep Field (CDF) North, but 
not the significantly larger correlation length observed in the CDF South. 

Subject headings: quasars: general - quasars: clustering - galaxies: evolution 



L INTRODUCTION 

Despite numerous observational efforts, quasar cluster- 
ing and its dependency on luminosity remains controversial. 
Early studies suggesting that clustering decreases with red- 
shift (lovino & Shaver 1988; Croom & Shanks 1996) are 
opposed by more recent observations, which suggest a more 
complicated clustering history that gradually increases with 
redshift (Kundic 1997; La Franca et al. 1998; Porcani et 
al. 2004; Croom et al. 2005). In particular, z < 2 studies 
with the Sloan Digital Sky Survey (SDSS, Myers et al. 2006, 
2007a, 2007b) and 2QZ (Croom et al. 2005, Porcani & Nor- 
berg 2006) have uncovered weak evidence for clustering evo- 
lution, while above z = 2 the SDSS indicates strong evolution 
in clustering (Shen et al. 2007) and 2QZ shows a somewhat 
weaker increase (Croom et al. 2005). Over all these redshifts, 
the luminosity dependence of clustering is weak (Porciani et 
al. 2004; Adelberger & Steidel 2005; Croom et al. 2005; My- 
ers 2006, 2007a), which is usually interpreted as being prob- 
lematic for quasar models in which AGN luminosity is corre- 
lated with proxies for halo mass. 

While questions of obscuration and completeness surround 
optical selection techniques, hard X-ray observations are 
largely unaffected by obscuration, thus making them per- 
haps the best candidate for identifying AGNs (Mushotsky 
2004). At redshifts z < 1 considerable effort has been put 
into measuring X-ray selected AGNs in both hard (2-10 keV) 
and soft (0.5-2 keV) bands (e.g. MuUis et al. 2004; GiUi et 
al. 2005; Basilakos et al. 2004, 2005; Yang et al. 2006; Miyaji 
et al. 2007). While optical z < 1 surveys tend to produce cor- 
relations lengths between 5-6 h~^ Mpc (e.g. Porciani & Nord- 
berg 2006), X-ray selected catalogs at these redshifts give 
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correlation lengths in the range 7-8 h ' Mpc (e.g. Mullis et 
al. 2004). However, the X-ray selected AGN from the Chan- 
dra Deep Field-North and South (CDF-N and CDF-S), exhibit 
significant variances in both the correlation length and power- 
law slope fits. Thus, Gilli et al. (2005), using the full catalog 
from 0.5-8 keV, found a correlation length of 5.5 ±0.6 h~^ 
Mpc for the CDF-N, in agreement with optical surveys, and 
10.3 ± 1.7 h'^ Mpc for the CDF-S, which is clearly in dis- 
agreement. This is surprising given that the log A^- log 5 of 
the two fields agree well. A recent re-analysis of this data, 
binned by luminosity and separating into hard and soft classi- 
fications, has yielded essentially the same overall result (Plio- 
nis et al. 2008), but elucidated that the two fields have much 
more consistent clustering behavior when binned by luminos- 
ity. Previous analyses had suggested the difference in cluster- 
ing could be attributed solely to sample variance due to the 
lack of large superclusters in the CDF-S (Gilli et al. 2003). 

While AGN/quasar feedback has become theoretically fa- 
vored as a necessary component of galaxy evolution (e.g. 
Scannapieco & Oh 2004, hereafter SO04; Granato et al. 2004; 
Croton et al. 2006), direct observational evidence of this idea 
remains somewhat weak. The paucity of high-redshift X-ray 
objects also makes studying clustering evolution of X-ray cat- 
alogs difficult. However, Francke et al. (2008) have cross- 
correlated AGN from the Extended CDF-S and luminous blue 
galaxy sources at z ~ 3 identified in the MUSYC survey (Ga- 
wiser et al. 2006). Their results indicate that the AGN tar- 
geted in the survey are more clustered than star-forming lumi- 
nous blue galaxies, a result consistent with the idea that typi- 
cal AGNs tend to sit in more massive halos than the average 
galaxy population. Clearly there remain many open questions 
about the clustering properties of optical and X-ray selected 
AGN, and the ongoing CDF controversy suggests that we can 
learn a great deal by comparing to predictions of clustering 
from simulations. 

From a theoretical perspective, early models of quasar for- 
mation associated quasars with galaxy mergers and assumed 
a close relationship between black hole mass and luminos- 
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ity (e.g. Kauffmann & Haehnelt 2000; Wyithe & Loeb 2002, 
2003). In this case, the black hole mass was calculated ei- 
ther by using the Mbh-ct relationship (Ferrarese & Mer- 
rit 2000) or associating Mbh with the halo circular velocity 
(Merrit & Ferrarese 2001; Tremaine et al. 2002; Ferrarese 
2002). Such "light bulb" models successfully match the lu- 
minosity function of high-redshift quasars (e.g. Wyithe & 
Loeb 2003), but become progressively more inaccurate at low 
redshifts when feedback processes become important (e.g. 
SO04). Cosmological simulations using the Mbh-o- frame- 
work and incorporating feedback have managed to reproduce 
the turn-down in the quasar luminosity function with moder- 
ate success (Thacker et al. 2006, hereafter TSC06). However, 
there appear to be differences between the detailed behavior 
of gas in simulations versus semi-analytic models, which are 
primarily due to differences between shock-heating in a uni- 
form medium relative to an inhomogeneous one (e.g. Helley 
et al. 2003; Nagamine et al. 2005; Croton et al. 2006; Catta- 
neo et al. 2006; Cattaneo et al. 2007; Di Matteo et al. 2008). 

More recent models (e.g. Hopkins et al. 2005a, 2005b, 
2005c, 2006a, 2006b, 2007a), motivated by numerical mod- 
eling of black hole accretion during mergers (Di Matteo et 
al. 2005; Springel et al. 2005a, 2005b; Robertson et al. 2006a, 
2006b; Cox et al. 2006a, 2006b), suggest that quasar activ- 
ity is comparatively decoupled from galaxy mass. This pic- 
ture entails complex relationships between a distinct sequence 
of AGN evolutionary epochs and the feedback processes that 
regulate them (e.g. Hopkins et al. 2007a). The resulting be- 
havior is one in which the bright end of the luminosity func- 
tion corresponds to quasars radiating at close to their peak 
luminosities near the Eddington limit, while the faint end cor- 
responds to the same population radiating in the faint part of 
their light curve, at or below w 0.1 of the Eddington lumi- 
nosity (Hopkins et al. 2005b). As a result, clustering is only 
a weak function of luminosity (Lidz et al. 2006). While the 
exact dynamics of nuclear accretion flows are still beyond the 
resolution of simulations of colliding galaxies, and are still the 
subject of much active research and modeling (e.g. Proga et 
al. 2008), the overall phenomenology in this model of quasar 
activity is weU understood. Recent increases in computing ca- 
pacity have lead to simulations of this model in a cosmologi- 
cal environment and investigations of the impact of AGN on 
disk formation (e.g. Sijacki et al. 2007; Di Matteo et al. 2008; 
Okamoto et al. 2008). 

In our earlier work (TSC06) we incorporated the merger 
and feedback AGN model of SO04 into a large cosmolog- 
ical smooth particle hydrodynamic simulation. However, 
our analysis of the clustering properties of optically-selected 
AGNs was constrained to a single luminosity bin and did not 
consider redshift evolution in any significant depth. In view 
of several recent ground-based surveys, it is therefore timely 
to reanalyze our simulation to address both the luminosity de- 
pendence and redshift evolution of clustering. Furthermore, 
the advent of new X-ray selected AGN catalogs with optical 
follow-up also allows us to present an analysis of the cluster- 
ing of an X-ray selected catalog. 

We stress that the aim of this paper is not to encourage sup- 
port for one quasar model over another, but rather to examine 
whether the details of the faint part of the light curve are actu- 
ally needed to accurately predict currently observed clustering 
statistics. Since the SO04 model does not include a low lumi- 
nosity accretion period we can indirectly constrain the impor- 
tance of this epoch to quasar clustering behavior. This should 
not be interpreted as constraining whether such a period does 



actually occur, or for that matter, the relative length of such a 
period. 

The layout of the paper is as follows: in §2 we give a sum- 
mary of our simulations and overall method. In §3 we present 
a detailed analysis of clustering at z = 1 .5 , 2.0, and 3 .0 and the 
dependence of this clustering on luminosity. In §4 we exam- 
ine the clustering properties of X-ray selected AGN, again as 
a function of redshift and luminosity. We close with a brief 
discussion in §5. Throughout the paper we consider a pre- 
WMAP3 (Spergel etal. 2003) ACDM model with parameters 
h = 0.7, flo = 0.3, Qa = 0.7, = 0.046, an = 0.9, and n = 1, 
where h is the Hubble constant in units of 100 km s"' Mpc^' , 
rto, Qa, and fib are the total matter, vacuum, and baryonic 
densities in units of the critical density, is the variance of 
linear fluctuations on the 8/!"'Mpc scale, and n is the "tilt" 
of the primordial power spectrum. While we consider only a 
fixed (78 in our simulation, the overall impact of changing erg 
on fitted correlation functions (^(r) = (ro/r)^) is to change the 
correlation length ro. For correlation functions with 7 « 2, in- 
creasing CTg by a factor of / will increase the correlation length 
of unbiased tracers of mass by the same factor Throughout 
the paper the Eisenstein & Hu (1999) transfer function is used 
and we quote all distances in comoving coordinates. 

2. SIMULATION METHOD AND QUASAR MODELING 

We consider two simulations in this study. The first is 
a "fiducial" run containing star formation and our model of 
AGN outflows (TSC06) in a periodic cube 146 /i"' Mpc on 
a side, containing 2 x 640^ particles. With these choices the 
dark-matter particle mass is 1.9 x IO^Mq and the gas parti- 
cle mass is 2.7 x 10^ Mq. The second simulation, which we 
call the "comparison" run, uses 2 x 320^ particles in a pe- 
riodic cube of size 73 li~^ Mpc, which matches the particle 
mass in the fiducial run, and includes star formation but not 
AGN outflows (Scannapieco et al. 2008). Both simulations 
were conducted with a parallel OpenMP based implementa- 
tion of the "HYDRA" code (Thacker & Couchman 2006) that 
uses the Adaptive Particle-Particle, Particle-Mesh algorithm 
to ciilculate gravitational forces (Couchman 1991), and the 
smooth particle hydrodynamic method to calculate gas forces 
(Lucy 1977; Gingold & Monaghan 1977). Due to computa- 
tional cost as well as the limitations of our modeling, both 
simulations were halted at z = 1 .2. 

Our method, as outlined in TSC06, associates quasar-phase 
AGN with galaxy mergers, which are tracked within the sim- 
ulations by identifying gas groups and applying group num- 
ber labels to their particles. Mergers are groups for which 
at least 30% of the accreted mass does not come from a sin- 
gle massive progenitor, and we calculate the mass of the as- 
sociated black hole, Mbh, using the circular velocity of the 
new system, Vc, and the observed Mbh - Vc relation (Merrit & 
Ferrarese 2001; Tremaine et al. 2002; Ferrarese 2002). We 
note that observational evidence of the universality of this re- 
lationship is still restricted to low redshift systems, and there 
is modest evidence that the normalization changes somewhat 
with redshift (e. g. Woo et al. 2008) . Continuing, our modeling 
approach yields 



Mbh = 2.8 x 10^ ( 
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and G is the gravitational constant, Pv(z) is the virial density 
as a function of redshift, and is the impHed virial radius for 
a group of gas particles with mass 



li/3 



4/37rpv(z) 



(3) 



In keeping with the model outlined in Wyithe & Loeb 
(2002) and SO04, we assume that for each merger the accret- 
ing black hole shines at its Eddington luminosity (1.2 x 10^** 
ergs s"' Mq~') for a time taken to be a fixed fraction, 0.055, 
of the dynamical time of the system, Tagn = 0.055rv/vc = 
5.8 X 10"''O(z)"'''^//(z)"'. We have shown in eariier work 
(TSC06) that apart from a small discrepancy at the most lu- 
minous end of the luminosity function, these simple assump- 
tions lead to a model that reproduces the observed AGN lu- 
minosity function as well as the predictions of SO04. Fur- 
thermore, we were able to demonstrate that this discrepancy 
can be explained in terms of the relative efficiency of shock 
heating on substructure, and corrected for if necessary by 
post-processing our simulations, as discussed in further detail 
below. An initial analysis of clustering properties was also 
in close agreement with observations, particularly the small- 
scale, r <; Mpc, clustering of quasars (e.g. Hennawi et 
al. 2006). 

Perhaps the most problematic aspect of this model is that 
it makes no distinction between AGN formed by gas-rich 
"wet" mergers versus those formed by gas-poor "dry" merg- 
ers. As we implement star formation on the basis of a simple 
merger model and ignore the quiescent mode of star forma- 
tion, it is difficult for us to make this distinction with any con- 
fidence. However, while there is significant evidence that at 
low redshifts the quiescent mode of star formation dominates 
(Noeske et al. 2007), at higher redshifts there is good reason 
to believe that mergers are necessary to fuel observed high 
star formation rates {e.g. Erb 2008). We can also appeal to the 
fact that while at z = dry mergers are important, they will 
be less so at z = 1 .2. For example, the semi-analytic estimates 
presented in Hopkins et al. (2007), in particular their Figure 
5, show that that the ratio of gas-rich to gas poor 4 x lO'^M© 
mergers, ranges from 10:1 at z = 2, to 2:1 at z = 1, indicating 
that gas-poor mergers are relatively unimportant before our 
final redshift. 

In the calculation of wind velocity we assume that a fixed 
fraction = 0.05 of the bolometric energy of each AGN is 
channeled into a kinetic outflow, while the remainder is emit- 
ted as light. While there is much debate about variability of 
this value on a system-by-system basis, our choice is consis- 
tent with other literature estimates (e.g. Furlanetto & Loeb 
2001; Nath & Roychowdhury 2002), as well as observations 
(Chartas et al. 2007). If we restrict ourselves to the con- 
sideration of systems with large bulges, then the resulting 
level of kinetic energy input is considerably greater than that 
from supernovae and stellar winds (e.g. Kravtsov & Yepes 
2000; Tozzi et aZ. 2000; Brighenti & Mathews 2001; Babul 
et al. 2002; Tornatore et al. 2004). Using the Eddington lu- 
minosity, associated dynamical time and the wind efficiency, 
each AGN outflow is launched with a wind energy of 
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Since there are considerable uncertainties about the precise 
geometry of AGN outflows, we have chosen to use a spherical 
shell to represent the outflow. Even strongly bipolar systems 



will tend to release an ellipsoidal cocoon of gas (e.g. Begel- 
man & Cioffi 1989; Yamada et al. 1999), so this approxima- 
tion is reasonable. We thus model each expanding outflow 
as a spherical shell at a radius 2ry which is created by rear- 
ranging the gas between and Ir^ which lies below a density 
threshold of 2.5py. In Figure [T]we present a plot of local gas 
density in a 12 Mpc region with the position of the virial 
and launching radii indicated. The radial velocity of the shell 
v, is set by ensuring that the sum of the thermal and kinetic 
wind energies is equal to /i^ - £^giav where figiav is the poten- 
tial energy change required to move the particles to 2ry. The 
post-shock temperature of the wind is given by 

This model produces a level of preheating in galaxy clusters 
and groups that is in good agreement with observations as dis- 
cussed in TSC06, to which the reader is referred for further 
details about our simulations. 

3. OPTICALLY-SELECTED AGN 

3.1. Optical Quasar Luminosity Function Revisited 

As a test of our overall approach, our first step is to re- 
peat the optical luminosity function analysis of TSC06, but 
for both simulations and for ranges of redshifts binned so as to 
make them most useful for comparisons with clustering mea- 
surements. As previously, we construct the luminosity func- 
tion by binning in luminosity and redshift. We calculate the 
number of quasars in each bin times the total time these ob- 
jects are shining, and divide by the time interval, the width of 
the bin, and the volume of the simulation. That is for a given 
redshift bin / and a given luminosity bin j the luminosity func- 
tion is simply 

1 NT^ 
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where the sum is over the lifetimes of all quasars with red- 
shifts and luminosities associated with the /, /' bin, which 
spans a time interval Af,- and a range of luminosities ALbj. 
The resulting luminosity functions for the fiducial AGN- 
feedback run and the comparison run are shown in Figure 
|2l in which the error bars are 1-sigma estimates, computed 
as A^ij = ±(l+NijT^^^], where Nij is the number 

of quasars contributing to bin i,j. As discussed in TSC06, 
our fiducial model shows a clear turn-down in the number 
of Lb > 10'^ quasars at z < 2, which parallels the observa- 
tional trend, but still overestimates the number of luminous 
and low-redshift quasars due to numerical effects. Likewise, 
the luminosity function in the no-feedback simulation con- 
tinues to rise at low redshift, increasing along with the halo 
merger rate as discussed in Wyithe & Loeb (2003). Thus, this 
simple model fails to reproduce the drop in the number den- 
sity of z < 2 quasars as discussed in SO04 and Scannapieco et 
al. (2008). 

Finally, we include a luminosity function calculated to 
match the shock behavior in the semi-analytic SO04 model. 
This was achieved by removing neighboring objects from the 
simulation that are found inside a shock radius calculated us- 
ing the SO04 model (see TSC06 for an extended discussion of 
this analysis). After we apply this algorithm, the simulation 
results much more closely match the observations. 

3.2. Dependence of the Correlation Function on Redshift and 

Luminosity 



Fig. 1 . — Schematic representation of tlie impact of outflows on the local density of gas. The region shown is 12 h ' Mpc across (comoving), at a redshift of 
z=1.59. The the virial radius, r,., of a system with a baryonic mass of 6 X IO'^Mq is represented by the inner circle, the outer circle corresponds to 2rv, denoting 
the launching point of an outflow for this system. A number of outflow events with different characteristic radii are visible within this small volume. 

Having outlined the successes and limitations of our model 
in reproducing the observed number density of quasars, we 
next move on to a detailed study of their spatial distribution. 
Here our primary tool is the real-space auto-correlation func- 
tion, calculated as 



^qq(r,Z,L)+l 



DD(r,z,L) 
RR{r,z,L) ' 
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where DD{r,z,L) is the number of pairs at a given comov- 
ing distance within a given redshift bin and with a luminosity 
within a given interval, and RR(r,z,L) is the average number 
of such pairs that would be found at this separation in a ran- 
dom distribution. Here we have correlated all quasars within 
each redshift bin, regardless of whether the two objects are 
shining simultaneously. While this vastly improves the sta- 
tistical signal, the use of a relatively large redshift window 
places a lower limit on the spatial scales that we can study, 
because the peculiar motions can shift the positions of the 
quasars during the finite time window associated with each 
bin. For our choices of redshift intervals, and estimating typi- 
cal peculiar velocities of quasars at w 300 km/s, this places 
a lower limit of 0.5 Mpc. Note that intrinsic velocity 
dispersion is estimated from the properties of the halos in 
which the majority of quasars are contained, and is somewhat 
smaller than observed pairwise dispersions (e.g. da Angela et 
al. 2005), which include both intrinsic and observational er- 



rors. 

The finite volume of our simulation places an upper limit 
on the distance we can study of approximately 1/5 the box 
size (Scoccimarro 1998; Szapudi et al. 1999), which corre- 
sponds to 30 Mpc in the AGN feedback run and 15 
Mpc in the comparison run. To be especially conservative 
we use a cutoff radius of lO/i"' Mpc for most of our results, 
which allows a direct comparison of the fiducial and com- 
parison runs. For the fiducial run alone we also examine the 
impact of changing to a 25/1"' Mpc cutoff. In each bin the 
error bars have be computed using a simple 1-sigma Poisson 
estimate of ACqq(r,z,L) = ^qq(r,z,L)[l ± (1 +DD(r,z,L)r^^'-] 

With these limitations in mind, in Figure [3] we plot the 
correlation function of simulated quasars, dividing our sam- 
ple into three luminosity bins from Lb = 10" - IO'^Lq.b, 
Lb= IO'^-IO'^Lq^b, and Lb = 10"- 10''^Lq,b. Focusing first 
on the AGN feedback run, the most striking feature of this 
plot is the relative lack of clustering at large separations and 
low redshifts, which occurs even though our model assumes 
accretion at the Eddington rate for all active black holes. Per- 
haps contrary to initial expectations, this weak dependence 
blankets the range of redshifts and separations that are best 
constrained observationally, suggesting that complex accre- 
tion histories may not play a key role in explaining current 
optical measurements. 
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Fig. 2. — Evolution of the B-band quasar luminosity function. The simulation results are given by the solid circles, while the dotted line is the simple estimate 
from the analytic model of Wyithe & Loeb (2003). From left to right the columns give results at redshifts of 1.2 — 1.75, 1.75-2.25, and 2.25-4.0. From top 
to bottom, the rows show results from the fiducial run, the fiducial ran with additional suppression imposed (by removing neighboring systems below a heating 
threshold), and the comparison ran. In all panels en'or bars are 1-sigma Poisson estimates. The observational data are taken from Croom et al. (2004, crosses) 
Richards et al. (2005, open triangles), Richards et al. (2006, open squares). Wolf et al. (2003, open circles), and Siana et al. (2008, filled squares). 



Note, however, that very different behavior occurs both 
at the smallest separations and at the highest redshifts, with 
^qq showing a strong luminosity dependence in both these 
regimes. Each of these enhancements is likely to be caused 
by different processes. At small separations, the strong depen- 
dence is likely to be a manifestation of so-called "one-halo ef- 
fects" (e.g. Beriind & Weinberg 2002; Bullock efaZ. 2002; van 
den Bosch et al. 2003; Magliocchetti & Porciani 2003) which 
take place when gravitationally -bound objects orbit each other 
within the same potential, adding significantly to the correla- 
tion function at distances smaller than the virial radius. This 
small-scale upturn, which has been confirmed observationally 
for SDSS quasars (Hennawi et al. 2006; Serber et al. 2006; 
Myers et al. 2007b), occurs at the radius corresponding to 
the maximum apocenter of such gravitationally-bound pairs, 
which in turn corresponds to the virial radius of the halos in 
which they are contained. As the most luminous AGN live 
in the deepest gravitational potential wells in our model, this 



means that small-scale clustering is naturally enhanced for 
these objects, causing a break in ^qq(r) that occurs at larger 
radii for more luminous objects. 

On the other hand, the luminosity dependence seen at ^ 2 

Mpc /;-' in the z = 3 bin is on such large scales that it can not 
be due to this effect. Instead this enhancement is likely to be a 
result of "geometrical bias," which is caused by the statistics 
of peaks within a Gaussian random field {e.g. Kaiser 1984; 
Bardeen et al. 1986; Mo & White 1996; Porciani et al. 1998). 
Note that on these scales the increase in Cqq('") is independent 
of distance, further pointing to this origin. 

Similar trends are apparent in the comparison run. Again 
there is little luminosity dependence at z < 2 and r > 1 /i"' 
Mpc. Also as in the feedback case, strong luminosity depen- 
dence is detected in the two regimes that are least constrained 
observationally: small-scale < 2 Mpc clustering, which 
is likely to be dominated by one-halo effects, and larger-scale 
high-redshift clustering, which is likely to be dominated by 
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Fig. 3. — Luminosity and redshift dependence of the quasar auto-correlation function, ^qq(r). Tlie top row sliows tlie results our AGN-feedback simulations, as 
calculated by partitioning quasars into bins with Lg = lO" - IO'^Lq.s, Lb = lO'" — IO'^Lq g, and Lg = lO'^^ - 1O''*L0 g. In the second row we show the ratio of 
5qq(r) for Lb = lO'^ - IO'^'L© quasars over 5qq(r) for Lg = lO'^ - IO'^Lq quasar (squares) and the ratio for Lb = lO" — IO'-Lq quasars over Lg = lO'^ - lO'^L©, 
quasars again from the AGN-feedback simulation. In the third row we show 5qq(r) from our no feedback comparison simulation, with symbols as above, and the 
ratios of the correlation functions from this ran are given the bottom row. As in Figure 2, from left to right each column shows the results from z= 1.2— 1.75, 
z= 1.75 — 2.25, and z = 2.25 -4.0, and all error bars are 1-sigma Poisson estimates. 



geometric bias. 

To quantify our results further we computed the bias of 
quasars as function of L and z as 



bHL,z) = 



(8) 



where ^qqifkjLjZ) is the quasar auto-correlation function in a 
radial bin k as a function of luminosity and redshift, ^DM{f,z) 
is the linear dark matter correlation function extrapolated to 
the redshift bin of interest, w{rk,L,z) is a weighting function 
that counts the number of pairs contributing to each value of 
^qqirk^L^z), and we average over the interval from r= 1.0 to 
10 /i"' Mpc. Error bars are again computed from Poisson esti- 
mates. A selected list of the computed bias values is given in 
Table[T]and the full data-set is plotted in Figure|4] in which we 
have also compiled results from current surveys, extrapolating 
to the B-band with a spectral slope of a^, = -0.5, (Wyithe & 
Loeb 2005; TSC06). Note however that these surveys do not 
necessarily estimate bias at exactly the same range of separa- 
tions as we have used. In particular, this plot includes points 
from Porciani et al. (2004) and Porciani & Norberg (2006), 
measured at sa 2-20 comoving Mpc, Adelberger & Stei- 



del (2005), measured at < 30 comoving h ' Mpc, Croom et 
al. (2005), measured at « 1-20 comoving /i"' Mpc, and My- 
ers et al. (2006; 2007a), measured from « 1 - 100 Mpc. 
Finally, for comparison purposes, we have also added a simple 
analytic estimate of the bias expected in the no feedback case, 
in which black hole mass can be directly related to the halo 
velocity dispersion and hence to the halo mass as in Wyithe & 
Loeb (2002). In this case. 



bmz)-- 



/2(l-c)_ 



i2c 



,i2c 



+ b(l-c)(l-c/2) 



(9) 

where a = 0.707, b = 0.5, c = 0.6 , 6o,, = 1.69, j/' = 
a^/^(5o,c£'(z)~'cr"'(Mhaio), D(z) is the hnear growth factor, 
o'(Mhaio) is the z = variance on the halo mass scale, Mhaio, 
corresponding to a given quasar luminosity (Sheth et al. 2001 ; 
see also Mo & White 1996; Jing 1999; Scannapieco & 
Barkana 2002), in the case in which gas accretion and dark 
matter collapse occur simultaneously, maintaining the cosmo- 
logical ratio at all times. 

Expressing our > 1 /i"' Mpc results as a bias allows for 
easy quantification of the trends seen in Figure [3] as well as 
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Fig. 4. — Top: Bias of quasars as a function of redshift and B-band luminosity. As in Figure 2, simulation results are given by the solid points with 1-sigma 
Poisson error bars, and the dashed Hne is the simple Sheth et al. (2001) estimate of the bias as described in the text. The open squares corresponds to the 
observational data points, which are taken from Porciani et al. (2004) and Porciani & Norberg (2006), open triangles; Croom et al. (2005), crosses; Adelberger 
& Steidel (2005), stars; and Myers et al. (2006; 2007a), squares. From top to bottom, the rows correspond to the fiducial run, the fiducial run with additional 
suppression, and the comparison run. Columns correspond to redshift bins as in Figures 2 and 3. 



comparisons with observations. Focusing first on the AGN 
feedback resuhs, we find that bias increases by no more than 
w 50% over the range of luminosities and redshifts probed 
by current surveys. The observational data indicate that cur- 
rent clustering bias measurements (Porcani et al. 2004; Adel- 
berger & Steidel 2005; Croom et al. 2005; Porcani & Norberg 
2006; Myers et al. 2006, 2007a) do not provide any significant 
statistical constraints above \ogyQ{LB / Lq) = 12.5. Indeed, a 
sample large enough to detect luminosity dependence of bias 
with A/? ~ 1 at a 3ct confidence level, given current detection 
limits (such as 2QZ), would require an all-sky measurement 
(Porcani & Norberg 2006). At z = 3, there is a suggestion in 
our simulations that bias changes significantly with luminos- 
ity above logiLs) = 12, but this regime is poorly-constrained 
observationally. 

While this mismatch between the observed range of red- 
shifts and luminosities and the regime in which we expect 
strong luminosity dependence appears initially to be some- 
what of a conspiracy, it can be understood naturally as a con- 
sequence of our feedback modeling. As discussed in detail 



in Scannapieco et al. (2005) and TSC06, AGN outflows act 
to impose a maximum halo mass, above which gas is unable 
to cool efficiently, suppressing further generations of galaxies 
and quasars. Furthermore, as radiative cooling is proportional 
to the square of the gas density, cooling is much more efficient 
in dense, high-redshift structures than it is at lower redshifts. 
This means that the "quenching threshold" i.e. the mass at 
which AGN is shut down (Faber et al. 2007) should decrease 
with time, with the strongest AGN quenching galaxy forma- 
tion even at high redshifts, but smaller smaller AGN quench- 
ing galaxy formation only at low redshifts. At the same time 
the hierarchical nature of dark-matter driven gravitational col- 
lapse means that the nonlinear mass scale increases with time, 
as ever-larger structures collapse and virialize. 

This means that AGN are naturally divided into two 
regimes. At high redshift, the characteristic luminosity of 
active black holes brightens along with the nonlinear mass 
scale, while at low redshift, their characteristic luminosity 
fades along with the quenching threshold. The z « 2 peak 
of AGN activity then marks a distinct transition between hier- 
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TABLE 1 

Selected correlation lengths and biases for the fiducial and comparison runs as a function of redshift, luminosity, and 
selection band. optically-selected values are given first. a key of 'f' corresponds to the fiducial run, while 'c to the 

comparison run. BOLOMETRIC luminosities associated WITH EACH BIN ARE GIVEN, ALONG WITH THE ASSOCIATED B-BAND LUMINOSITY (FOR 
THE OPTICAL CATALOG) AND HARD X-RAY LUMINOSITY (FOR THE X-RAY CATALOG). CORRELATION LENGTHS ARE QUOTED TO 2 SIGNIFICANT 
FIGURES WITHOUT ERRORS SINCE STATISTICAL ERRORS WILL BE SIGNIFICANTLY SMALLER THAN SYSTEMATIC ERRORS FROM THE BINNING 

PROCEDURE. 



Run 


z 


log(LBo;/L0) 


log(LB/L0) 


Mb 


Lx 
(erg s"') 


cut 

Mpc) 


(r' Mpc) 


b 


F 


1.5 


12.7 


11.6 


-24.3 




10 


5.2 


1.9 


F 


1.5 


12.7 


11.6 


-24.3 




25 


4.4 


1.6 


c 


1.5 


12.7 


11.6 


-24.3 




10 


4.1 


1.5 


F 


2.0 


12.7 


11.6 


-24.3 




10 


5.3 


2.3 


F 


2.0 


12.7 


11.6 


-24.3 




25 


4.6 


2.0 


c 


2.0 


12.7 


11.6 


-24.3 




10 


3.8 


1.6 


F 


3.0 


12.7 


11.6 


-24.3 




10 


5.4 


3.1 


F 


3.0 


12.7 


11.6 


-24.3 




25 


5.3 


3.0 


c 


3.0 


12.7 


11.6 


-24.3 




10 


4.3 


2.5 


F 


1.5 


13.2 


12.1 


-25.6 




10 


6.4 


2.4 


F 


1.5 


13.2 


12.1 


-25.6 




25 


6.3 


2.2 


c 


1.5 


13.2 


12.1 


-25.6 




10 


4.8 


1.8 


F 


2.0 


13.2 


12.1 


-25.6 




10 


7.2 


3.0 


F 


2.0 


13.2 


12.1 


-25.6 




25 


6.2 


2.6 


c 


2.0 


13.2 


12.1 


-25.6 




10 


5.6 


2.4 


F 


3.0 


13.2 


12.1 


-25.6 




10 


7.1 


4.0 


F 


3.0 


13.2 


12.1 


-25.6 




25 


7.4 


4.0 


c 


3.0 


13.2 


12.1 


-25.6 




10 


4.6 


2.7 


F 


1.5 


13.4 


12.4 


-26.2 




10 


7.8 


3.1 


F 


1.5 


13.4 


12.4 


-26.2 




25 


5.7 


2.1 


c 


1.5 


13.4 


12.4 


-26.2 




10 


6.3 


2.3 


F 


2.0 


13.4 


12.4 


-26.2 




10 


8.6 


3.5 


F 


2.0 


13.4 


12.4 


-26.2 




25 


6.9 


2.9 


c 


2.0 


13.4 


12.4 


-26.2 




10 


6.7 


3.0 


F 


3.0 


13.4 


12.4 


-26.2 




10 


9.1 


5.1 


F 


3.0 


13.4 


12.4 


-26.2 




25 


8.3 


4.4 


C 


3.0 


13.4 


12.4 


-26.2 




10 


7.5 


4.1 


F 


1.75 


10.0 






3.2 X lO''^ 


10 


3.9 


1.6 


C 


1.75 


10.0 






3.2 X 10^*2 


10 


3.1 


1.2 


F 


1.75 


11.3 






3.2 X 10^*3 


10 


3.6 


1.5 


C 


1.75 


11.3 






3.2 X 10"*^ 


10 


3.2 


1.3 


F 


3.0 


11.3 






3.2 X 10"*^ 


10 


3.3 


2.0 


C 


3.0 


11.3 






3.2 X 10"*^ 


10 


3.2 


1.8 


F 


1.75 


12.6 






3.2 X lO"*"* 


10 


5.2 


2.1 


C 


1.75 


12.6 






3.2 X lO"*"* 


10 


4.5 


1.7 


F 


3.0 


12.6 






3.2 X lO"*"* 


10 


5.5 


3.1 


C 


3.0 


12.6 






3.2 X 10'*4 


10 


4.5 


2.5 



archical and anti-hierarchical formation, which occurs when 
the quenching threshold drops below the nonlinear mass scale. 
Thus the majority of AGN formed at redshifts below the peak 
of AGN activity, meaning those that are easiest to observe and 
quantify, are naturally found in halos with masses well below 
the nonlinear mass scale. As can been seen from eq. (|9]l, these 
low masses are very weakly biased, as is a strong function 
of i^' when ly' ^l, but almost a constant when i^' < I. 

Another important feature of our feedback model is that it 
produces halo biases that are systematically offset from the 



no-feedback case. This is because feedback acts to slow ac- 
cretion even before the quenching threshold is passed, mean- 
ing that each halo hosts a somewhat less massive black hole 
than it would have in the absence of feedback. Thus, for a 
fixed luminosity, each AGN is shifted to a somewhat more 
massive, and hence more clustered, dark matter halo, and the 
typical increase in bias over the no-feedback case is about 
30%. At the faint end this corresponds to an increase over the 
analytic estimate of about a factor of « 2. The quasars points 
from our simulation are also offset from the simple Sheth et 
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Fig. 5. — Top: Correlation length as a function of luminosity and redshift. 
The closed triangles, circles, and squares correspond to redshift bins z= 1.5, 
2.0, and 3.0 respectively, with 1-sigma Poisson error bars. Open squares 
con'espond to the 7 = 1.8 fits of Porciani & Norberg (2006) to the 2QZ survey 
at redshifts z = 0.93, 1.19, 1.41, 1.60, 1.79, and 1.98, with the mean Mb being 
taken from Croom et al. (2005). Open triangles correspond to the Shen et 
al. (2007) SDSS data points for 2.9 < z < 3.5 and z > 3.5, but are fitted for 
7 = 2.0 rather than 1.8 with i-band to B-band conversion taken from Hao et 
a!. (2005). 

al. (2005) estimates, again because gas accretion lags behind 
the dark matter collapse in the simulation. 

It is important to note that post-processing our results to 
correct for the handful of very luminous, low-redshift AGN 
that result from inefficiency in shock heating in the simula- 
tion has very little effect on any of these trends, despite it 
leading to an almost perfect match of the luminosity function. 
As shown in the center row, removing these objects only im- 
pacts the Lb ^ lO'^L© measurements in the lowest redshift 
bins, primarily increasing the already large error bars by fur- 
ther lowering the number density of these objects. 

In the comparison run, on the other hand, the overall bias at 
each mass scale is somewhat lower than in either of the other 
two cases. Again this is because while gas accretion and cool- 
ing still take time, this time is much less than the 2-5 Gyr 
Hubble times at these epochs, meaning that the relationship 
between black hole mass and halo mass is more in line with 
that expected purely from the dark matter distribution. Note 
that even in this run however, there is very little evolution in 
clustering over the observed range of luminosities. This is be- 
cause even though no feedback is included, leading to a fair 
number of large and biased low-redshift AGN in the simula- 
tion, the lack of low redshift quasars in the data in the obser- 
vations means these objects simply do not exist in nature, and 
thus can not be compared to our predictions. However, even 
in the absence of feedback, the significant cooling time asso- 
ciated with large objects means that gas accretion trails dark 
matter collapse. Thus means that quasars are found in higher- 
mass halos, and hence are more clustered, than one would 
expect in a simple model in which the gas accretion moves 
forward in lock-step with dark-matter collapse. 

A second way to quantify our results is by using the corre- 



lation length, ro, the scale at which = 1. If cx r this 
occurs at = £,qq{r)r'^ for all choices of r. Thus we can com- 
pute ro by averaging this quantity over all radial bins from 
1 to 10 /r' Mpc 



ro(L,z) = 



1/7 



(10) 



where we choose to set 7 = 1.8 and ^qq is the quasar auto- 
correlation function as a function of luminosity L and redshift 
z. In Figure |5] we compare the results of this analysis with 
correlation length measurements from the 2QZ (range 1 < z < 
2) and SDSS (range 3 < z < 5) surveys. A selected list of the 
computed correlation lengths is also given in Table[T] 

The fiducial model, shown in the top panel, indicates that 
below Lb ~ lO'^L© there is little evolution in the correlation 
length either as a function of redshift or luminosity. The cal- 
culated correlation lengths are also in good agreement with 
the Croom et al. (2005) results. Above Lb ~ lO'^L©, on the 
other hand, the variability of the correlation with luminosity 
is more evident, however even this seems to fall slightly short 
of the very large correlation length seen in the SDSS sam- 
ple (although our quoted Poisson errors are smaller than the 
systematic errors from the binning procedure) . We note that 
our our highest redshift bin (z = 2.25-4.0) has a mean red- 
shift less than that of the z > 3.5 SDSS bin, but the lower 
2.9 < z < 3.5 SDSS bin is comparable to ourz = 3 predictions. 
The primary difference between results probably stems from 
the fact that the SDSS measurement considers pairs with sep- 
arations from 4/1"' Mpc < r < 150/!~'Mpc, which is a vastly 
larger range than can be probed with current simulations that 
retain high resolution in individual galaxies. In fact, it could 
be reasonably argued that a box of 1000 Mpc is neces- 
sary to study correlations on this scale. Additionally, consid- 
ering correlations to such large radii is fraught with potential 
difficulties, since for redshifts z < 1 the galaxy-galaxy corre- 
lation function is expected to steepen for separations larger 
than 60/!"' Mpc (see Springel, Frenk & White 2006). Were 
AGN/quasars to more closely trace the underlying dark mat- 
ter auto-correlation function than the normal galaxy popula- 
tion at high redshift then the effective power law for the auto- 
correlation function would be expected to change by almost 
50% as one moves from 5 Mpc to 20 /i"' Mpc (Porciani & 
Norberg 2006). Further, fitting steeper power laws will inher- 
ently tend to produce larger correlation lengths. 

The steepening of the galaxy-galaxy correlation function 
is a product of the underlying bias of the galaxy population 
and the transition to the non-linear regime in the dark matter 
power spectrum, which occurs at A: > ~ 0. 1 /j Mpc"' at low 
redshift. This scale evolves comparatively slowly until z ~ 1, 
above which it begins to recede quickly, and by z — 3 we find 
k„i ~ 0.3 /; Mpc"' . Semi-analytic models of the galaxy-galaxy 
correlation function are able to roughly preserve the location 
of this steepening point (see, e.g., Springel et al. 2006) but do 
so through a rapidly increasing bias with redshift. If the bias 
of the AGN population does not increase sufficiently quickly 
with redshift then the clustering statistics will directly mea- 
sure the underlying evolution in the non-linear scale. It is 
worth noting that current surveys at z ~ 3 definitely straddle 
the turnover in the dark matter correlation function. Thus the 
enhanced correlation seen in the SDSS might be related to a 
change in the radial position of the steepening of the correla- 
tion function with redshift. The precise details are dependent 
upon the redshift evolution of the bias of the AGN population 
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and, as we have indicated, our simulation box is too small to 
make any firm statements. 

However, to study the dependence of our results on the 
much smaller distances we can probe, we recalculated eq. 
(flOb using a range of separations from r^^ = 1 /;"' Mpc to 25 
/z"' Mpc. These values, shown in the second panel of Fig- 
ure |5] demonstrate the correlation is decreased systematically 
when one includes more information from large separations. 
As discussed above in relation to Figure |3] the origin of these 
differences is most likely to be the excess contribution at small 
separations from the one-halo term, which becomes less im- 
portant as moves from the shorter 10 Mpc cutoff to the 
longer 25 /;"' Mpc cutoff. Even this modest change in the 
outer cutoff can change the correlation lengths by as much as 
40% (specifically in the z=1.5,Lb = 12.4 bin) although the 
mean change is close to 15%. 

Finally, in the lower panel, we show the results from our 
comparison simulation, which again displays similar trends as 
in the AGN feedback run, but with a lower level of clustering. 
For all models and separations it is difficult to draw conclu- 
sions about the redshift evolution of ry. Indeed, our results 
support the idea of using correlation function evolution mod- 
els that are roughly constant in comoving coordinates since 
we see little evolution in our comparison models (for exam- 
ple) at log,o(LB/Lboi) < 12.5. 

4. X-RAY SELECTED AGN 

X-ray selection is widely believed to be an unbiased method 
for selecting AGN candidates, which is largely free from the 
obscuration and incompleteness issues that affect optical cat- 
alogs (e.g. Yang et al. 2006). While the exact nature of the 
optical versus X-ray light curves is the subject of debate (as 
summarized in Hopkins et al. 2008, and references therein), 
we examine the X-ray clustering of our catalog based on the 
same assumptions as our optical catalog. In this case the only 
differences between X-ray and optical selection come from 
bolometric correction factors, as even obscuration of optical 
systems should not impact the overall correlation unless ob- 
scuration is somehow a function of position. Additionally, so 
long as the average lifetime of the X-ray bright period is not 
longer than w 1 Gyr, and hence peculiar velocities do not lead 
to systems moving an appreciable distance, our lifetime as- 
sumptions should not have a significant impact on clustering. 
Of course this is not true for the luminosity function, which is 
very sensitive to changes in lifetimes, and we have previously 
shown (TSC06) that the hard X-ray luminosity function cal- 
culated from our model reproduces the observations of Ueda 
et al. (2003). 

To calculate the ratio of the intrinsic X-ray luminosity, 
Lx, to the bolometric luminosity, LboI, within our simulation 
we use the following two-polynomial fits from Marconi et 
al. (2004), 

log[LBol/Lz(2-10keV)] = 1.54+0.24£+0.012£2+0.0015/:^ 

log[LBoi/ix(0.5-2keV)] = 1.65+0.22£+0.012£2+0.0015£^ 

(11) 

where C = log(LBoi)- 12, and LboI is given in Lq. While a 
large number of estimates for bias and correlation lengths ex- 
ist for z <; 1 (e.g. MulUs et al. 2004; Gilli et al. 2005; Basi- 
lakos et al. 2004, 2005; Yang et al. 2006; Miyaji et al. 2007), 
the current observational data is too sparse at redshifts z > 2 
to provide reliable statistics on redshift evolution. How- 
ever, cross-correlating luminous blue galaxies and AGN al- 
lows a calculation of the bias of the AGN population at z ~ 3 
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Fig. 6. — Bias of X-ray selected AGNs as a function of redshift and X- 
ray luminosity. The green triangles are from the comparison run, while the 
blue circles are from the fiducial feedback model, both with 1-sigma Pois- 
son error bars. The open square data point on the z= 1-V5 plot is from Yang 
et al. (2006), and con'esponds to the bias for their z = 1.5 — 3.0 bin (no vari- 
ance for the mean of the luminosity bin is given). The open triangle is from 
Francke et al. (2008) and gives their AGN bias calculated from the cross- 
correlation function of AGN and luminous blue galaxies at z ~ 3. The error 
bars encapsulate their range in luminosity and the mean is likely lightward of 
the central value. While our fiducial run is in good agreement with the Yang 
et al. (2006) results, the z — 3 result is clearly lower than Francke et al. (2008) 
data. 

(Francke et al. 2008). We therefore calculate both the bias 
and correlation functions of our simulation, but split the cata- 
log into only two redshift bins: z = 1.2-2.0, which we label 
as z = 1.75, and z = 2.0-4.0, which we label as z = 3.0. Within 
these bins we then calculate correlation functions and bias us- 
ing the procedures outlined in section 3.2. 

In Figure|6]we plot the bias of our catalog for both hard and 
soft bands in our two redshift bins. As for the optical cata- 
log, we find that the simulation with quasar feedback has a 
higher bias than the comparison run. Again, the primary ori- 
gin of this difference is that feedback forces AGN of a given 
luminosity into more massive, and thus more biased halos. 
The fiducial run hard X-ray data in the z = 1 .75 bin are a 
good fit to the Yang et al. (2006) bias estimate from the CDF- 
N and CLASXS fields, which suggests that including feed- 
back is necessary to bias halos sufficiently to match observa- 
tions. Furthermore, comparison to the z = 3 bias estimates of 
Francke et al. (2008), calculated using the cross-correlation 
of AGN and luminous blue galaxies, shows that our numbers 
are low relative to these observations, and that AGN feedback 
is absolutely necessary to match these data. We achieve a 
marginal agreement if we take the low value of their error and 
also the higher end of the luminosity range, which is plausible 
since we plot the center of the range and the the mean of their 
bin is likely rightward of the central value. It is also worth 
noting that these results show that in our model the luminos- 
ity dependence begins to become more noticeable above 10"^^ 
erg s"'. 

Before examining our results further we mention that since 
the luminosities of the X-ray data we are considering are con- 
siderably lower than the equivalent optically-selected catalog, 
the observations we consider may well be probing different 
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Fig. 7. — Correlation function of X-ray selected AGNs. We have divided 
the simulation into two redshift ranges, z = 2.25-4.0, which we label z = 3, 
andz =1.1 -2.25 which we label z = 1.75. We then bin into 3 decades of 
Lx, ranging from lO'" erg s"', to 10'*'' erg s"'. The ratio of the 3.2 X 10'''' 
erg s"' to 3.2 X 10'"' erg s"' X -ray correlation functions increases by a factor 
of three in the one-halo regime, with a similar rise being observed in the 
optically selected catalogs. All error bars are 1-sigma Poisson estimates. 

AGN fueling mechanisms as compared to our major merger 
model. In particular Hopkins & Hernquist (2006) have sug- 
gested that the limiting upper luminosity for where secular 
(i.e. Seyfert) effects or minor mergers become important is 
around Mb — -22. This corresponds to an X-ray luminosity 
of 10'*'* erg s"', which is roughly in the middle of our con- 
sidered X-ray luminosity range. However, since we compare 
directly to observations with mean luminosities above 10"''* 
erg s"' we can be reasonably confident that major mergers are 
the dominant physical process in these systems. 

Both the fiducial and comparison runs correspond to the 
same bias at the faint end, ignoring the faintest bins at z = 3 
which are at the limits of our resolution. However, the overall 
sensitivity to luminosity is higher for the fiducial run when ex- 
amined over the entire range Lx = [10'*'', lO"*^'^] erg s"' , show- 
ing an increase in the bias of roughly a factor of 4 in the hard 
and soft X-ray bands at z = 1 .75, and a factor of close to 6 in 
the z = 3.0 bin. Hints of this dependence are observed in the 
optical catalog over the range logiLs) = [1 1 .8, 12.6], although 
it is difficult to determine visually since the brighter bins have 
large error bars while the faint-end cutoff at Lb = lO^ '^L© 
does not probe as low in luminosity as the X-ray catalog. In 
the X-ray data, the ratio of the hard X-ray bias of the fidu- 
cial run to the comparison run is w 1 .6 for Lx = 10*^'^ erg s"' 
at both z= 1.75 and z = 3.0 and the soft X-ray numbers are 
similar. The Lx = lO"*^^ erg s~' bin also shows redshift evolu- 
tion, with its bias decreasing by a factor of two from z = 3 to 
z= 1.75. 

In Figure [T] we plot the correlation function of our X-ray 
selected AGN. While our estimates seem to be in good agree- 
ment with the data (the little studied z = 3 bin aside), the 
correlation lengths in this plot appear to be smaller than the 
observed data. Table [T] quantifies the correlation lengths ex- 
tracted from the X-ray data using eq. (fTOl l. as well as selected 
correlation lengths and biases from throughout this paper Ob- 
servationally, Yang et al. (2006) give a combined CDF-N, 



CLASXS X-ray correlation length of ro = 6.l!y;^/i"' Mpc, al- 
beit with a shallow 7 = 1.47 slope. Matching this correlation 
length at z = 1 .75 with our luminosity binned data requires a 
mean luminosity greater than 3.2 x 10"*"* erg s"*. However, 
just like the optical data these surveys have a much higher 
cutoff radius than our simulation. For example the CLASXS 
field discussed in Yang et al. (2006) considers pairs with sep- 
arations up to 200 /;"' Mpc, well beyond beyond the radius at 
which the down-turn occurs in the galaxy-galaxy correlation 
function. 

Recently, using luminosity binning, Plionis et al. (2008) 
have suggested that the wide variation in the correlation 
lengths of the CDF-S and CDF-N can be reconciled. They 
find evidence of strong evolution in clustering as a function 
of luminosity, with the correlation length in the hard X-ray 
band increasing from w 6/i~' Mpc to w 18/;"' Mpc, with a 0.7 
dex increase in luminosity. A comparison to our X-ray cor- 
relation lengths in Table [T] shows that even our fiducial "light 
bulb" model cannot produce this level of luminosity depen- 
dence, nor can we produce the same underlying clustering. 
Thus although we do find a somewhat stronger trend for ro to 
increase with Lx than Lb, our results are at odds with Plionis et 
al. (2008) It is also worth noting that their quoted correlation 
length for the highest flux in the soft band of w 30 Mpc 
is considerably larger than that quoted for the IRAC Shallow 
Cluster Survey (Brodwin et al. 2007), ro = 19.14!^;^^/;-' Mpc 
at z = 0.97. This implies that the comparatively low lumi- 
nosity AGN (mean Lx ~ lO^-' erg s"') they sample are more 
strongly clustered at high redshift than z = 1 galaxy clusters. 
Ultimately more clustering data is needed to help understand 
the high redshift clustering of X-ray selected AGN, and we 
eagerly anticipate future all-sky surveys. 

5. DISCUSSION 

We have presented an analysis of our simulated 
quasar/AGN catalog, focusing on the dependence of 
real-space clustering on redshift, luminosity, and selection. 
Our model does not follow the detailed accretion history 
onto the central supermassive black hole (e.g. Hopkins et 
al. 2005a, 2005b, 2005c, 2006a, 2006b, 2007a), but rather 
takes a simple one-to-one correspondence between black 
hole mass and luminosity. Nonetheless we capture much of 
the essential physics in AGN formation and feedback. As 
a consistency check on our earlier work, we showed that a 
model that does not include feedback precisely follows the 
predicted luminosity function of the Wyithe & Loeb (2003) 
model, and thus fails badly at low redshift by overpredicting 
the observed numbers counts. 

On the other hand, the qualitative luminosity function be- 
havior is reproduced (TSC06) when feedback is included. 
Our clustering results are also in broad agreement with ob- 
served data, the main difference being somewhat less evolu- 
tion with redshift than observed and somewhat smaller cor- 
relations lengths, although our bias values are in quantitative 
agreement within the observational errors. Furthermore for 
our "light bulb" model, the dependence of clustering with lu- 
minosity is weak at the luminosities probed by current sur- 
veys. Although the underlying relationship between quasar 
luminosity and black hole mass is likely to more complex than 
the simplified model assumed in this study (e.g. Ganguly et 
al. 2007), modeling these complexities dose not appear to be 
necessary to understanding current clustering measurements. 

While the assumptions used in our calculations limit our 
analysis to scales above w 0.5 Mpc we are still able to 
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clearly observe one-halo effects, which occur within w 2 Mpc 
for luminous AGN. Significantly, the luminosity depen- 
dence is more visible in this part of the correlation function as 
the two systems are embedded in a more highly biased halo. 
For example, the ratio of the 3.2 x 10'*^' erg s"' to 3.2 x 10'*'* 
erg s"' X-ray correlation functions increases by a factor of 
three in the one-halo regime. Deep quasar pair data would 
thus be extremely useful in helping to determine luminosity 
dependence in more detail. 

However, we reemphasize that a straightforward compari- 
son of current observations to our results, or those of any sim- 
ulations, is not possible. While larger samples have made ob- 
servational studies more robust, there still remain differences 
in fitted scales for correlation functions and the assumed slope 
of power-law fits. This is particularly important in the context 
of calculating correlation lengths, as departures from power- 
law fits occur both on small scales and large scales. On small 
scales the one-halo term produces a steepening in the effec- 
tive index, while on large scales the transition from the non- 
linear to linear regime in the dark-matter power spectrum also 
produces a steepening of the effective index. We also note 
that redshift-space distortions (e.g. Croom et al. 2005; da An- 
gela et al. 2005, 2008), can also produce departures from pure 
power-law behavior Ultimately, the true power-law slope ob- 
served for the AGN/quasar population will depend on the un- 
derlying bias. 

Given these facts and our modest outer radius of 25 
Mpc, the fits we derive should be treated with due caution. For 
example, even for modest changes in the outer cutoff from 25 
Mpc to 10 /i"' Mpc, our correlation lengths can increase 
by as much as 40%, although the mean change is close to 
15%. The increase is directly associated with the one-halo 
contribution being given more weight in the case with the 
shorted outer cutoff, although we emphasize that all fits be- 
low 10 Mpc are within the non-linear scale at the epochs 
we are considering. 



The redshift and luminosity dependence of large-scale clus- 
tering is a product of two competing effects: growth of the 
non-linear mass scale with time and a decrease in the mass 
scale of the quenching threshold that limits the supply of fuel 
to AGN. This quenching is primarily a function of the mean 
density of the gas, which controls its cooling rate. Thus at 
high redshifts, when radiative cooling is extremely efficient, 
feedback is weak, and the luminosity of black holes grows 
along with the nonlinear mass scale. However, once feed- 
back is able to heat gas to a cooling time longer than the Hub- 
ble time, the fuel supply for luminous AGN is quenched, and 
this quenching become more efficient as the universe expands. 
The net result is a peak in AGN activity at redshift z w 2: 
AGN formed after this redshift correspond to low-mass, low- 
bias halos and show a weak luminosity dependence, and AGN 
formed before this redshift correspond to a wide range of bi- 
ases and show an appreciable luminosity dependence. Thus 
the clustering behavior of AGN is a direct result both of the 
evolution of dark matter halos and the physics of AGN feed- 
back. 
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